Welcome to the third lesson of AeroPython! We created some very interesting potential flows in lessons 1 and 2, with our Source & Sink notebook, and our Source & Sink in a Freestream notebook.
Think about the Source & Sink again, and now imagine that you are looking at this flow pattern from very far away. The streamlines that are between the source and the sink will be very short, from this vantage point. And the other streamlines will start looking like two groups of circles, tangent at the origin. If you look from far enough away, the distance between source and sink approaches zero, and the pattern you see is called a doublet.
Let's see what this looks like. First, load our favorite libraries.
In [1]:
import math
import numpy
from matplotlib import pyplot
In the previous notebook, we saw that a source-sink pair in a uniform flow can be used to represent the streamlines around a particular shape, named a Rankine oval. In this notebook, we will turn that source-sink pair into a doublet.
First, consider a source of strength $\sigma$ at $\left(-\frac{l}{2},0\right)$ and a sink of opposite strength located at $\left(\frac{l}{2},0\right)$. Here is a sketch to help you visualize the situation:
The stream-function associated to the source-sink pair, evaluated at point $\text{P}\left(x,y\right)$, is
$$\psi\left(x,y\right) = \frac{\sigma}{2\pi}\left(\theta_1-\theta_2\right) = -\frac{\sigma}{2\pi}\Delta\theta$$Let the distance $l$ between the two singularities approach zero while the strength magnitude is increasing so that the product $\sigma l$ remains constant. In the limit, this flow pattern is a doublet and we define its strength by $\kappa = \sigma l$.
The stream-function of a doublet, evaluated at point $\text{P}\left(x,y\right)$, is given by
$$\psi\left(x,y\right) = \lim \limits_{l \to 0} \left(-\frac{\sigma}{2\pi}d\theta\right) \quad \text{and} \quad \sigma l = \text{constant}$$Considering the case where $d\theta$ is infinitesimal, we deduce from the figure above that
$$a = l\sin\theta$$$$b = r-l\cos\theta$$$$d\theta = \frac{a}{b} = \frac{l\sin\theta}{r-l\cos\theta}$$so the stream function becomes
$$\psi\left(r,\theta\right) = \lim \limits_{l \to 0} \left(-\frac{\sigma l}{2\pi}\frac{\sin\theta}{r-l\cos\theta}\right) \quad \text{and} \quad \sigma l = \text{constant}$$i.e.
$$\psi\left(r,\theta\right) = -\frac{\kappa}{2\pi}\frac{\sin\theta}{r}$$In Cartesian coordinates, a doublet located at the origin has the stream function
$$\psi\left(x,y\right) = -\frac{\kappa}{2\pi}\frac{y}{x^2+y^2}$$from which we can derive the velocity components
$$u\left(x,y\right) = \frac{\partial\psi}{\partial y} = -\frac{\kappa}{2\pi}\frac{x^2-y^2}{\left(x^2+y^2\right)^2}$$$$v\left(x,y\right) = -\frac{\partial\psi}{\partial x} = -\frac{\kappa}{2\pi}\frac{2xy}{\left(x^2+y^2\right)^2}$$Now we have done the math, it is time to code and visualize what the streamlines look like. We start by creating a mesh grid.
In [2]:
N = 50 # Number of points in each direction
x_start, x_end = -2.0, 2.0 # x-direction boundaries
y_start, y_end = -1.0, 1.0 # y-direction boundaries
x = numpy.linspace(x_start, x_end, N) # creates a 1D-array for x
y = numpy.linspace(y_start, y_end, N) # creates a 1D-array for y
X, Y = numpy.meshgrid(x, y) # generates a mesh grid
We consider a doublet of strength $\kappa=1.0$ located at the origin.
In [3]:
kappa = 1.0 # strength of the doublet
x_doublet, y_doublet = 0.0, 0.0 # location of the doublet
As seen in the previous notebook, we play smart by defining functions to calculate the stream function and the velocity components that could be re-used if we decide to insert more than one doublet in our domain.
In [4]:
def get_velocity_doublet(strength, xd, yd, X, Y):
"""Returns the velocity field generated by a doublet.
Arguments
---------
strength -- strength of the doublet.
xd, yd -- coordinates of the doublet.
X, Y -- mesh grid.
"""
u = - strength/(2*math.pi)*((X-xd)**2-(Y-yd)**2)/((X-xd)**2+(Y-yd)**2)**2
v = - strength/(2*math.pi)*2*(X-xd)*(Y-yd)/((X-xd)**2+(Y-yd)**2)**2
return u, v
def get_stream_function_doublet(strength, xd, yd, X, Y):
"""Returns the stream-function generated by a doublet.
Arguments
---------
strength -- strength of the doublet.
xd, yd -- coordinates of the doublet.
X, Y -- mesh grid.
"""
psi = - strength/(2*math.pi)*(Y-yd)/((X-xd)**2+(Y-yd)**2)
return psi
Once the functions have been defined, we call them using the parameters of the doublet: its strength kappa and its location x_doublet, y_doublet.
In [5]:
# computes the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)
# computes the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)
We are ready to do a nice visualization. To get the plots inside the notebook, we use the command %matplotlib inline :
In [6]:
%matplotlib inline
# plots the streamlines
size = 10
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u_doublet, v_doublet,
density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o');
Just like we imagined that the streamlines of a source-sink pair would look from very far away. What is this good for, you might ask? It does not look like any streamline pattern that has a ptactical use in aerodynamics. If that is what you think, you would be wrong!
A doublet alone does not give so much information about how it can be used to represent a practical flow pattern in aerodynamics. But let's use our superposition powers: our doublet in a uniform flow turns out to be a very interesting flow pattern. Let's first define a uniform horizontal flow.
In [7]:
u_inf = 1.0 # freestream speed
Remember from our previous lessons that the Cartesian velocity components of a uniform flow in the $x$-direction are given by $u=U_\infty$ and $v=0$. Integrating, we find the stream-function, $\psi = U_\infty y$.
So let's calculate velocities and stream function values for all points in our grid. And as we now know, we can calculate them all together with one line of code per array.
In [8]:
u_freestream = u_inf * numpy.ones((N, N), dtype=float)
v_freestream = numpy.zeros((N, N), dtype=float)
psi_freestream = u_inf * Y
Below, the stream function of the flow created by superposition of a doublet in a free stream is obtained by simple addition. Like we did before in the Source & Sink in a Freestream notebook, we find the dividing streamline and plot it separately in red.
The plot shows that this pattern can represent the flow around a cylinder with center at the location of the doublet. All the streamlines remaining outside the cylinder originated from the uniform flow. All the streamlines inside the cylinder can be ignored and this area assumed to be a solid object. This will turn out to be more useful than you may think.
In [9]:
# superposition of the doublet on the freestream flow
u = u_freestream + u_doublet
v = v_freestream + v_doublet
psi = psi_freestream + psi_doublet
# plots the streamlines
size = 10
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.contour(X, Y, psi, levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')
# calculates the stagnation points
x_stagn1, y_stagn1 = +math.sqrt(kappa/(2*math.pi*u_inf)), 0
x_stagn2, y_stagn2 = -math.sqrt(kappa/(2*math.pi*u_inf)), 0
# adds the stagnation points to the figure
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2], color='g', s=80, marker='o');
What is the radius of the circular cylinder created when a doublet of strength $\kappa$ is added to a uniform flow $U_\infty$ in the $x$-direction?
You have the streamfunction of the doublet in cylindrical coordinates above. Add the streamfunction of the free stream in those coordinates, and study it. You will see that $\psi=0$ at $r=a$ for all values of $\theta$. The line $\psi=0$ representes the circular cylinder of radius $a$. Now write the velocity components in cylindrical coordinates, find the speed of the flow at the surface. What does this tell you?
A very useful measurement of a flow around a body is the coefficient of pressure $C_p$. To evaluate the pressure coefficient, we apply Bernoulli's equation for ideal flow, which says that along a streamline we can apply the following between two points:
$$p_\infty + \frac{1}{2}\rho U_\infty^2 = p + \frac{1}{2}\rho U^2$$We define the pressure coefficient as the ratio between the pressure difference with the free stream, and the dynamic pressure:
$$C_p = \frac{p-p_\infty}{\frac{1}{2}\rho U_\infty^2}$$i.e.,
$$C_p = 1 - \left(\frac{U}{U_\infty}\right)^2$$In an incompressible flow, $C_p=1$ at a stagnation point. Let's plot the pressure coefficient in the whole domain.
In [10]:
# computes the pressure coefficient field
cp = 1.0 - (u**2+v**2)/u_inf**2
# plots the pressure coefficient field
size = 10
pyplot.figure(figsize=(1.1*size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
contf = pyplot.contourf(X, Y, cp, levels=numpy.linspace(-2.0, 1.0, 100), extend='both')
cbar = pyplot.colorbar(contf)
cbar.set_label('$C_p$', fontsize=16)
cbar.set_ticks([-2.0, -1.0, 0.0, 1.0])
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')
pyplot.contour(X,Y,psi, levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2], color='g', s=80, marker='o');
Show that the pressure coefficient distribution on the surface of the circular cylinder is given by
$$C_p = 1-4\sin^2\theta$$and plot the coefficient of pressure versus the angle.
Don't you find it a bit fishy that the pressure coefficient (and the surface distribution of pressure) is symmetric about the vertical axis?
That means that the pressure in the front of the cylinder is the same as the pressure in the back of the cylinder. In turn, this means that the horizontal components of forces are zero.
We know that, even at very low Reynolds number (creeping flow), there is in fact a drag force. The theory is unable to reflect that experimentally observed fact! This discrepancy is known as d'Alembert's paradox.
Here's how creeping flow around a cylinder really looks like:
In [11]:
from IPython.display import YouTubeVideo
YouTubeVideo('rbMx2NMqyuI')
Out[11]:
If you look carefully, there is a slight asymmetry in the flow pattern. Can you explain it? What is the consequence of that?
Here's a famous visualization of actual flow around a cylinder at a Reynolds number of 1.54. This image was obtained by S. Taneda and it appears in the "Album of Fluid Motion", by Milton Van Dyke. A treasure of a book.
In [12]:
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom.css', 'r').read()
return HTML(styles)
css_styling()
Out[12]: